library(ggplot2)

######################################
## Plot histograms of first saccades after onset of linguistic event of interest
## by region created by jdegen on 06/18/2013 Uses a reduced version of the
## dataset reported in Degen & Tanenhaus (under review). In particular, the
## garden-path and number term conditions are excluded. The dataset is
## pre-downsampled to 20ms.
# set your path here
load("processed_data/dt.RData")

# choose early or late subset
early = subset(dt, POD == "early")
late = subset(dt, POD == "late")

######################################
## exclusion criteria (taken from Salverda, Kleinschmidt, and Tanenhaus, under
## review): a) no saccade was generated after the onset of the target word (2.3%
## in SKT, 5% in DT); b) first saccade following the onset of the target word
## was directed to a region that the participant was already fixating (23.9% of
## remaining data in SKT, 38.9% in DT); c) first saccade following onset of
## target word was not directed to either target or competitor (50.6% in SKT,
## 26.3% in DT) in the DT early dataset, 44% of trials are _not_ excluded
## (compared to 36.8% in SKT)
######################################

######################################
# Step 1. Get all the first saccades
######################################

delay = subset(early, Time_rel_stim_Qonset > 0)
delay = subset(delay, Right_Type == "SACC")
tmp = delay[,c("TrialDataID","DataID")]
tmp = droplevels(tmp)
# get unique sample ID of first saccade
t = tapply(as.character(tmp$DataID),tmp$TrialDataID,FUN="[",1)
nrow(t)
row.names(delay) = delay$DataID
delay = delay[t,]


######################################
# Step 2. Exclude trials based on SKT criteria
######################################

# % trials excluded because there were no new saccades: 5%
1-(nrow(t)/length(unique(early$TrialDataID))) 

## if you want to remove trials where subject did NOT click on target
#table(delay$ResponseType)
#delay=subset(delay,delay$ResponseType=="t")

## add information about look on the sample directly prior to the first saccade (prevSampleRegion) and at word onset (FixWordOnset)
early$prevSampleRegion = as.factor(c(NA,as.character(early[1:(nrow(early)-1),]$rp_RegionType)))
row.names(early) = early$DataID
wordonset = subset(early, Time_rel_stim_Qonset == 0)
wordonset = droplevels(wordonset)
delay$prevSampleRegion = early[as.character(delay$DataID),]$prevSampleRegion
wordonset = subset(wordonset, TrialDataID %in% delay$TrialDataID)
delay$FixWordOnset = wordonset$rp_RegionType
head(delay,20)
nrow(delay) # 2318
save(delay,file="rdata/saccades_early.RData")
# subejct distribution
table(delay$FixWordOnset)
table(delay$SubjectName)
table(delay$SubjectName,delay$FixWordOnset)
t = as.data.frame(table(delay$SubjectName))
colnames(t) = c("Subject","Freq")
# histogram of how many subjects contributed how many trials
ggplot(t,aes(x=Freq)) +
  geom_histogram()


## exclude trials where first saccade is to a region the subject is already fixating
out.already=delay[as.character(delay$rp_RegionType)==as.character(delay$FixWordOnset),]
delay=delay[as.character(delay$rp_RegionType)!=as.character(delay$FixWordOnset),]
# excludes 38.9% of trials!!
nrow(out.already)/(nrow(delay) + nrow(out.already))

## exclude trials where first saccade was not to target or competitor
out.other = subset(delay,rp_RegionType %in% c("center","ci","ti","none"))
delay = subset(delay,rp_RegionType %in% c("t","c"))
delay = droplevels(delay)
# another 26.3% gone
nrow(out.other)/(nrow(delay) + nrow(out.other))
# leaves 1043 trials
nrow(delay)

########################################## 
# Step 3. Generate saccade histograms for each condition
########################################## 

# size of bins in histograms
HISTO_BINWIDTH = 40  # try other bindwidths too

delay$Region = delay$rp_RegionType
delay$Time = delay$Time_rel_stim_Qonset

ggplot(delay, aes(x = Time, fill = Region)) +
  geom_histogram(position = "identity", alpha = .5, binwidth = HISTO_BINWIDTH) +
  xlim(0, 800) +
  facet_grid(Quantifier~Numbers) +
  xlab("Time since quantifier onset (ms)") + ylab("Number of saccades") +
  theme_bw()


######################################
## Step 4. Moving window analysis to determine first time window in which probability of saccade to target is different from probability of saccade to competitor
######################################

chi_critical = 3.841 # two-tailed (?)
chi_marginal = 2.706

 